
  # Please always run steps 1-3 below before executing the code for specific Tables or Figures
  # Please set the working directory before starting
  # setwd()


  # 1. Upload libraries ##################################################################################    

  library(lfe)
  library(ggplot2)
  library(gridExtra)
  library(dplyr)
  library(data.table)
  library(boot) 
 
  ########################################################################################################   
  

  # 2. Functions #########################################################################################    

  # Cleans the data for each relevant variable
  .Clean.Data <- function(data,Variable,Pov,Effect1,New){
    
    if(New == 1){data <- subset(data,  Year > 2004  & !is.na(PABV))}
    if(New == 0){data <- subset(data,  Year <= 2004  )}
    data$Var <- data[,match(Variable,names(data))] 
    data <- subset(data, !is.na(Var))
    data <- subset(data, Poverty > Pov)

    
    ### for the variables that describe challenger's profiles
    if(Variable == "chl.entr.client" | Variable == "chl.entr.left" | Variable == "chl.entr.educ")
    {data <- subset(data, Inc.1 == 0 & Term.2 == 1) 
    }else{
      
      if(Variable == "chl.succ.client" | Variable == "chl.succ.left" | Variable == "chl.succ.educ")
      {data <- subset(data, Inc.1 == 0 & Term.2 == 1 & Rank <= 2) 
      }else{
        
        if(Variable == "pct" | Variable == "MV"  
           | Variable == "Ncand" | Variable == "Turnout"){data <- subset(data, Inc.1 == 1 & MV < 100)
        }else{
          
          data <- subset(data, Term.1 > 0) 
          data$d <- duplicated(data[,c(match("City",names(data)),match("Year",names(data)))])
          data <- subset(data, d == 'FALSE') 
          data$d <- NULL
        }
      }
    }
    
 
    data$Eff1 <- data[,match(Effect1,names(data))]
    data$Eff2 <- data$Year
    
    data$d <- duplicated(data[c("City","Year","Party")])
    data <- subset(data, d == 'FALSE')
    
    
    return(data)
  } 

  # Cleans the data for each relevant variable
  .Formulas <- function(IV,controls,Variable,New){
      fe <- " | Eff1+Eff2 "
      base1 <- "Var ~   D + Pop2 + IDH2 + D_Pop2 + D_IDH2  "
      base2 <- "Var ~   Pop2 + IDH2 + D_Pop2 + D_IDH2 "
      cont2 <- ""
      cont00 <- "+ LON*LAT + Area"
      cont <- " + LON*LAT  + Area + ESF + lGDPpc + Men + BF2_pt + OldCCT "
      if(Variable == "pct" ){cont2 <- " + pct.p"}
      if(Variable == "POOR_4" ){cont2 <- " +POOR00_4"}
      if(controls == 1){cost_f <- paste(cont,cont2,sep="")}else{cost_f <- ""}
      if(IV == 0) {formula <- as.formula(paste(base1,cost_f,fe,"|0|City",sep=""))}
      if(IV == 1) {formula <- as.formula(paste(base2,cost_f,fe,"| (BF2 ~ D) |City",sep=""))}
      if(New == 0 & (Variable == "pct") ){cont2 <- " + pct.p"}
      if(New == 0 & controls == 1){formula <- as.formula(paste(base1,cont00,cont2,fe,"|0|City",sep=""))}
      
      return(formula)
  }

  # Runs the regression for n of the 19 points in the treatment frontier
  .Run.coeff <- function(full.data,formula,sdPOP,sdIDH,kernel,bd,pop.cut,reg.beg,reg.end,pretre,IV){
    #### get the individual regressions
    tot <- 19
    
    x <- round((4/pi)**0.5,2) 
    bd <- bd*x
    mn1 <- mean(log(full.data$Pop)) 
    mn2 <- mean(full.data$IDH)                                      # averages for density
    sd1 <- sd(log(full.data$Pop)) 
    sd2 <- sd(full.data$IDH) ;                                        # sds for density
    cor12 <- cor(log(full.data$Pop),full.data$IDH)                                                    # correl for density
    vec <- matrix(0,tot,10) 
    
    med <- (tot + 1)/2 
    div <- 1 
    
    ii <- 1
   
    for (ii in reg.beg:reg.end){
      if(ii <= med){ pt2 <- 0.7  
                     pt1 <- (pop.cut-25)+ii*2.5*div}else
        { pt1 <- pop.cut 
          pt2 <- 0.7-0.0125*div*(ii-med) }       # set cut-off points       
      
      data <- data.frame(full.data )
      data$Pop2 <- data$Pop - pt1 
      data$IDH2 <- data$IDH - pt2                   # SET new variables
      data$Pop3 <- data$Pop2/sdPOP 
      data$IDH3 <- data$IDH2/sdIDH
      v1 <- as.numeric(bd[ii,1])
      v2 <- as.numeric(bd[ii,2])
 
      
      data <- subset(data, abs(Pop2) <= v1*sdPOP &  abs(IDH2) <= v2*sdIDH )         #apply bandwidth
      data$h <- (data$Pop3**2+data$IDH3**2)**0.5 
      sin <- as.vector( abs(data$Pop3)/data$h )
      cos <- as.vector( abs(data$IDH3)/data$h )
      data$r <- (v1*v2)/((( (v1*sin)**2 )+(  (v2*cos)**2  ))**0.5)
      data <- subset(data, h <= r)
      
      ## apply kernel
      if(kernel > 0){
        if(kernel == 1){data$wght1 <- v1 - abs(data$Pop3)}
        if(kernel == 2){data$wght1 <- 0.75*(1 - ((data$Pop3/v1)**2))}
        if(kernel == 3){data$wght1 <- (1/((2*pi)**2))*exp(-0.5*((abs(data$Pop3/v1)**2)))}
        data$wght1[data$D == 1] <- data$wght1[data$D == 1] / sum(data$wght1[data$D == 1])
        data$wght1[data$D == 0] <- data$wght1[data$D == 0] / sum(data$wght1[data$D == 0])
        if(kernel == 1){data$wght2 <- v2 - abs(data$IDH3)}
        if(kernel == 2){data$wght2 <- 0.75*(1 - ((data$IDH3/v2)**2))}
        if(kernel == 3){data$wght2 <- (1/((2*pi)**2))*exp(-0.5*((abs(data$IDH3)/v2)**2))}
        data$wght2[data$D == 1] <- data$wght2[data$D == 1] / sum(data$wght2[data$D == 1])
        data$wght2[data$D == 0] <- data$wght2[data$D == 0] / sum(data$wght2[data$D == 0])
        data$wght <- data$wght2*data$wght1
        data$wght[data$D == 1] <- data$wght[data$D == 1] / sum(data$wght[data$D == 1])
        data$wght[data$D == 0] <- data$wght[data$D == 0] / sum(data$wght[data$D == 0])
      }else{data$wght <- 1/nrow(data)}
      

      
      data$D_Pop2 <- data$D*data$Pop2 
      data$D_IDH2 <- data$D*data$IDH2   # reposition forcing variables
      data <- data.frame(data)
      
      if(pretre==0){
        xx <- felm(formula, data = data, weights = data$wght)# run regression
        r <- data.frame(summary(xx)$coefficients)
        lab <- "D"
        if(IV==1) lab <- "`BF2(fit)`"
      }
      
      if(pretre==1){
        
        data <- subset(data, D == 0)
        xx <- felm(Var ~   Pop2 + IDH2 |0|0|City, data = data, weights = data$wght)                                                         # run regression
        r <- summary(xx)$coefficients
        lab <- "(Intercept)"
        
      }
      
      vec[ii,1] <- r[lab,1] 
      vec[ii,2] <- r[lab,2]  
      vec[ii,3] <- r[lab,3] 
      vec[ii,4] <- nrow(data) 
      vec[ii,5] <- .Dnorm(log(pt1),pt2,mn1,mn2,sd1,sd2,cor12)
      vec[ii,6] <- summary(xx)$r.squared
      vec[ii,7] <- v1/x 
      vec[ii,8] <- v2/x
      vec[ii,9] <- pt1 
      vec[ii,10] <- pt2
     
    }
    
    vec <- data.frame(vec) 
    vec$low <- vec$X1 - 1.65*vec$X2 
    vec$high <- vec$X1 + 1.65*vec$X2
    return(vec)
  }
  
  # Creates the average of the several estimators
  .Run.coeff2 <- function(vec,reg.beg,reg.end){
    output <- matrix(0,1,10)

    vecx <- vec[1:19,] 
    output[1,1] <- weighted.mean(vecx[reg.beg:reg.end,1],vecx[reg.beg:reg.end,5])
    output[1,8] <- weighted.mean(vecx[reg.beg:reg.end,4],vecx[reg.beg:reg.end,5]) 
    output[1,9] <- weighted.mean(vecx[reg.beg:reg.end,7],vecx[reg.beg:reg.end,5])
    output[1,10] <- weighted.mean(vecx[reg.beg:reg.end,8],vecx[reg.beg:reg.end,5])    
    
    return(output)
  }
  
  # For bootstrapping standard errors
  .Boot.coeff <- function(base, data, indices, formula, std1, std2, IV, 
                          bd, kernel, pop.cut, reg.beg, reg.end, pretre){
    
    library(lfe)
    x <- round((4/pi)**0.5,2) 
    bd <- bd*x
    cit <- data[indices,]
    full.data <- merge(cit,base, by="City")
    mn1 <- mean(log(full.data$Pop)) 
    mn2 <- mean(full.data$IDH)                                      # averages for density
    sd1 <- sd(log(full.data$Pop)) 
    sd2 <- sd(full.data$IDH)                                          # sds for density
    cor12 <- cor(log(full.data$Pop),full.data$IDH)                                                    # correl for density
    output <- matrix(0,1,1)
    
    vec <- matrix(0,19,3) 
    vec[,3] <- 1 
    
    for (ii in reg.beg:reg.end){
      if(ii <= 10){pt2 <- 0.7 ;pt1 <- (pop.cut-25)+ii*2.5}else{pt1 <- pop.cut 
      pt2 <- 0.7-0.0125*(ii-10) } 
      
      v1 <- as.numeric(bd[ii,1])
      v2 <- as.numeric(bd[ii,2])
      
      a <- (1/(2*pi*sd1*sd2*((1-cor12**2)**0.5)))
      b1 <- ((log(pt1) - mn1)**2)/(sd1**2)
      b2 <- ((pt2 - mn2)**2)/(sd2**2)
      b3 <- -((log(pt1) - mn1)*(pt2 - mn2)*2*cor12)/(sd1*sd2)
      b <- exp(-(1/(2*((1-cor12**2)**0.5)))*(b1+b2+b3))
      vec[ii,2] <- (a*b)                      # joint density at the point
      
      data <- full.data 
      data$Pop2 <- data$Pop - pt1 
      data$IDH2 <- data$IDH - pt2             # SET new variables
      data <- subset(data, abs(Pop2) <= v1*std1 &  abs(IDH2) <= v2*std2 )         # apply bandwidth
      data$Pop3 <- data$Pop2/std1 
      data$IDH3 <- data$IDH2/std2 
      data$h <- (data$Pop3**2+data$IDH3**2)**0.5 
      data$sin <- abs(data$Pop3)/data$h 
      data$cos <- abs(data$IDH3)/data$h 
      x <- (4/pi)**0.5
      data$r <- (v1*v2)/((((v1*data$sin)**2)+((v2*data$cos)**2))**0.5)
      data <- subset(data, h <= r)
      
      if(kernel >= 1){
        if(kernel == 1){data$wght1 <- v1 - abs(data$Pop3)}
        if(kernel == 2){data$wght1 <- 0.75*(1 - ((data$Pop3/v1)**2))}
        if(kernel == 3){data$wght1 <- (1/((2*pi)**2))*exp(-0.5*((abs(data$Pop3/v1)**2)))}
        data$wght1[data$D == 1] <- data$wght1[data$D == 1] / sum(data$wght1[data$D == 1])
        data$wght1[data$D == 0] <- data$wght1[data$D == 0] / sum(data$wght1[data$D == 0])
        if(kernel == 1){data$wght2 <- v2 - abs(data$IDH3)}
        if(kernel == 2){data$wght2 <- 0.75*(1 - ((data$IDH3/v2)**2))}
        if(kernel == 3){data$wght2 <- (1/((2*pi)**2))*exp(-0.5*((abs(data$IDH3)/v2)**2))}
        data$wght2[data$D == 1] <- data$wght2[data$D == 1] / sum(data$wght2[data$D == 1])
        data$wght2[data$D == 0] <- data$wght2[data$D == 0] / sum(data$wght2[data$D == 0])
        data$wght <- data$wght2*data$wght1
        data$wght[data$D == 1] <- data$wght[data$D == 1] / sum(data$wght[data$D == 1])
        data$wght[data$D == 0] <- data$wght[data$D == 0] / sum(data$wght[data$D == 0])
      }else{data$wght <- 1/nrow(data)}
      
      data$D_Pop2 <- data$D*data$Pop2 
      data$D_IDH2 <- data$D*data$IDH2
      
      if(pretre==0){
        xx <- felm(formula, data = data, weights = data$wght)
        lab <- "D"
        if(IV==1) {lab <- "`BF2(fit)`"}
        vec[ii,1] <- summary(xx)$coefficients[lab,1]
      }else{
        datax <- subset(data, D==0)
        form <- as.formula("Var ~ Pop2 + IDH2  | 0 | 0 | 0")
        xx <- felm(form, data = datax, weights = datax$wght)
        vec[ii,1] <- summary(xx)$coefficients["(Intercept)",1]
      }
      
    }
    vec[,2] <- vec[,2]*vec[,3]
    output[1,1] <- weighted.mean(vec[reg.beg:reg.end,1],vec[reg.beg:reg.end,2])  
    
    return(output)
  }
  
  # Calculates BCA-adjusted confidence Intervals
  .Coeff.intervals <- function(output,results){
    ctt <- 1  
    for (pp in 1:nrow(output)){                         # by confidence level
      if(R >= 4000){ci <- boot.ci(results, conf = c(0.99,0.95,0.90),type="bca",index=pp)$bca} else{ci <- matrix(0,3,5)}
      for (ppp in 1:3){          # by coefficient 
        if (ppp == 1){ct <- 2} ; if (ppp == 2){ct <- 4} ; if (ppp == 3){ct <- 6}
        output[ctt,ct] <- ci[ppp,4]   
        output[ctt,ct+1] <- ci[ppp,5]
        
      }
      ctt <- ctt + 1
    }
    return(output)
  }
    
  # Calculates density of sample around a point in the frontier
  .Dnorm <- function(x1, x2, mn1, mn2, sdd1, sdd2, cor12){
    a <- (1/(2*pi*sdd1*sdd2*((1-cor12**2)**0.5)))
    b1 <- ((x1 - mn1)**2)/(sdd1**2)
    b2 <- ((x2 - mn2)**2)/(sdd2**2)
    b3 <- -((x1 - mn1)*(x2 - mn2)*2*cor12)/(sdd1*sdd2)
    b <- exp(-(1/(2*((1-cor12**2)**0.5)))*(b1+b2+b3))
    return(a*b)
  }
  
  # Creates table for pasting in LaTEX
  .LatexTable <- function(data,rownames,spaces){
    cols <- ncol(data)*2+1
    res <- data.frame(matrix(0,nrow(data),cols))
    res[,] <- "&"
    fill <- seq(from=2,to=cols,by=2)
    res[,fill] <- data
    res <- cbind(rownames,res)
    lastone <- rep("\tabularnewline",nrow(res))
    lastone[spaces] <- "\tabularnewline[0.25cm]"
    lastone[nrow(res)] <- "\tabularnewline[0.25cm]"
    res[,ncol(res)] <- lastone
    names(res) <- NULL
    return(res)
  }  
 
  #######################################################################################################   
 
 
  # 3. Open Files #######################################################################################    
  
  ### For all Tables and Figures
      M.data <- fread("Main_data.csv")
      M.data <- data.frame(M.data)
      
    ### Set SDs for scaling
      uni <- unique(data.table(M.data), by="City")
      sdPOP <- round(sd(uni$Pop),0) 
      sdIDH <- round(sd(uni$IDH),2) 
      rm(uni)

    ### Set treatment points
      pop.cut <- 30
      hdi.cut <- 0.7
      M.data$D <- 0  
      M.data$D[M.data$Pop <= pop.cut & M.data$IDH <= hdi.cut] <- 1 

    ### For ALL
      band.matrix <- data.frame(fread("optimal_bands.csv"))
 
    ### For Figure II only
      dex <- fread("Coverage_data.csv")
    
  #######################################################################################################   
  
  
  
  # Table I #############################################################################################    
  
    reg.beg <- 9                              # First point of the 19 in the frontier to be estimated
    reg.end <- 14                             # Last point of the 19 in the frontier to be estimated

    vars <- c("lPABV","BF2")

    R <- 5000
    matr <- data.frame(matrix(0,4,6))
  
    i <- 2 ; j <- 1
    
    for (i in 1:length(vars)){
      for(j in 0:5){
  
        set.seed(55)
  
        # Set column specifications for Table
        if(j<4)           bd <- subset(band.matrix, Var == vars[i] & Kernel == 1)[,1:2] 
        if(j==0)          {pretre <- 1}else{pretre <- 0}
        if(j %in% c(2,3)) {ctr <- 0}else{ctr <- 1}
        if(j==4)          {bd <- matrix(.9,19,2)}
        if(j==5)          {bd <- matrix(.75,19,2)}
        
        ### Cleans Data and Run regressions   
        data <- .Clean.Data(M.data,Variable=vars[i],Pov=25,Effect1="State",New=1)
        data$x <- 1
        formula <- .Formulas(IV=0,controls=ctr,Variable=vars[i],New=1) 
        if(j==2) formula <- as.formula("Var ~ D + Pop2 + IDH2 + D_Pop2 + D_IDH2 |0|0|City")

        vec <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                                        kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                                         reg.end=reg.end, pretre=pretre, IV=0) 
        output <- .Run.coeff2(vec,reg.beg,reg.end) 
        cit2 <- data.frame(unique(data$City),1) 
        names(cit2) <- c("City","One") 
       
        ### Bootsrap SDs
        results <- boot(data=cit2, statistic=.Boot.coeff, R=R, base = data, std1 = sdPOP , 
                                        std2= sdIDH , kernel = 1, formula = formula,
                                        bd = bd, IV = 0, parallel = "snow", 
                                        ncpus = 8,reg.end = reg.end, reg.beg = reg.beg, 
                                        pop.cut = pop.cut, pretre=pretre)
        
        output <- .Coeff.intervals(output,results)  
        
        ### prepare output
        OUT <- data.frame(output) 
 
        names(OUT) <- c("COEFF","BCA99L","BCA99H","BCA95L","BCA95H","BCA90L","BCA90H","Obs","Bpop","Bhdi") #; OUT <- cbind(aux,OUT) 
        OUT$Stars <- ""
        if(j > 0){
        OUT$Stars[(OUT$BCA90L*OUT$BCA90H) > 0] <- "*" ; 
        OUT$Stars[(OUT$BCA95L*OUT$BCA95H) > 0] <- "**" 
        OUT$Stars[(OUT$BCA99L*OUT$BCA99H) > 0] <- "***"
        }
        
        coeffi <- OUT$COEFF[1]
        bcalow <- OUT$BCA90L[1]
        bcahig <- OUT$BCA90H[1]
        
        if(i==1 & j==0){
          coeffi <- exp(coeffi)/1000000
          bcalow <- exp(bcalow)/1000000
          bcahig <- exp(bcahig)/1000000
        }
        
        matr[1,j+1] <- paste(format(round(coeffi, 3),nsmall=3),OUT$Stars[1],sep="")
        matr[2,j+1] <- paste("[",format(round(bcalow, 2),nsmall=2),",",
                           format(round(bcahig, 2),nsmall=2),"]",sep="")
        matr[3,j+1] <- paste("{",format(round(OUT$Obs[1], 0),nsmall=0),"}",sep="")
        matr[4,j+1] <- paste("(",format(round(OUT$Bpop[1], 2),nsmall=2),",",
                           format(round(OUT$Bhdi[1], 2),nsmall=2),")",sep="")
        print(j)
    
      }
      
      if(i==1) {matrx <- matr}else{matrx <- rbind(matrx,matr)}
      
    }
    
    ccc <- c("Health funds, R$bn","CI","Obs. per bin",
             "Bandwidths","CCT cov, p.p.","CI ","Obs. per bin ","Bandwidths ")
    table <- .LatexTable(matrx,ccc,c(4,8))

    print(table, row.names=FALSE)
  
  #######################################################################################################   
    
  # Table II ############################################################################################    
    
    reg.beg <- 9                              # First point of the 19 in the frontier to be estimated
    reg.end <- 14                             # Last point of the 19 in the frontier to be estimated
    ctr <- 1                                  # Include controls
   
    vars <- c("pct","MV","Ncand","POOR_4",
              "chl.entr.educ","chl.entr.client",
              "chl.succ.educ","chl.succ.client")
    
    R <- 5000
    matr <- data.frame(matrix(0,2,5))
    
    for (i in 1:length(vars)){
      for(j in 0:3){
        
        set.seed(55)
        
        # Set column specifications for Table
        bd <- subset(band.matrix, Var == vars[i] & Kernel == 1)[,1:2] 
    
        if(j==0)          {pretre <- 1}else{pretre <- 0}
        if(j==2)          {NN <- 0}else{NN <- 1}
        if(j==3)          {IV <- 1}else{IV <- 0}
        
        ### Cleans Data and Run regressions   
        data <- .Clean.Data(M.data,Variable=vars[i],Pov=25,Effect1="State",New=NN) 
        
        formula <- .Formulas(IV=IV,controls=ctr,Variable=vars[i],New=NN) 
        
        vec <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                          kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                          reg.end=reg.end, pretre=pretre,IV=IV) 
        output <- .Run.coeff2(vec,reg.beg,reg.end) 
        cit2 <- data.frame(unique(data$City),1) 
        names(cit2) <- c("City","One") 
        
        ### Bootsrap SDs
        results <- boot(data=cit2, statistic=.Boot.coeff, R=R, base = data, 
                        std1 = sdPOP , 
                        std2= sdIDH , kernel = 1, formula = formula,
                        bd = bd, IV = IV, parallel = "snow", 
                        ncpus = 8,reg.end = reg.end, reg.beg = reg.beg, 
                        pop.cut = pop.cut, pretre=pretre)
        
        output <- .Coeff.intervals(output,results)  
        
        ### prepare output
        OUT <- data.frame(output) 
        
        names(OUT) <- c("COEFF","BCA99L","BCA99H","BCA95L","BCA95H",
                        "BCA90L","BCA90H","Obs","Bpop","Bhdi") 
        OUT$Stars <- ""
        if(j > 0){
          OUT$Stars[(OUT$BCA90L*OUT$BCA90H) > 0] <- "*" ; 
          OUT$Stars[(OUT$BCA95L*OUT$BCA95H) > 0] <- "**" 
          OUT$Stars[(OUT$BCA99L*OUT$BCA99H) > 0] <- "***"
        }
        
        matr[1,j+1] <- paste(format(round(OUT$COEFF[1], 3),nsmall=3),OUT$Stars[1],sep="")
        matr[2,j+1] <- paste("[",format(round(OUT$BCA90L[1], 2),nsmall=2),",",
                             format(round(OUT$BCA90H[1], 2),nsmall=2),"]",sep="")
        
        if(j==3) { matr[1,5] <- paste("(",format(round(OUT$Bpop[1], 2),nsmall=2),",",
                                             format(round(OUT$Bhdi[1], 2),nsmall=2),")",sep="")
        matr[2,5] <- paste("{",format(round(OUT$Obs[1], 0),nsmall=0),"}",sep="") }
        
        print(j)
        
      }
      
      if(i==1) {matrx <- matr}else{matrx <- rbind(matrx,matr)}
      
    }
    
    ccc <-   c("Incumbent's vote share","(%)","Margin of victory","(p.p.)",
                      "Candidates","(number)","Pro-poor spending,","(% share)",
                           "Challenger's entry ","(share with HS)","Challenger's entry",
               "(share Clientelistic)","Challenger is top 2 ","(share with HS) ",
               "Challenger is top 2","(share Clientelistic) ")
    
    table <- .LatexTable(matrx,ccc,c(16))
    
    print(table, row.names=FALSE)
    
  #######################################################################################################   
    
  # Table III ###########################################################################################    
    

    parties_client <- c(11,12,14,15,22,25) 
    parties_nonclient <- c(13,40,65,23,45) 
    parties_left <- c(13,12,16,40,65,23,50,29,21,18,33) 
   
    Pov <- 25                                 # Minimum poverty level (%)
    reg.beg <- 9                              # First point of the 19 in the frontier to be estimated
    reg.end <- 14                             # Last point of the 19 in the frontier to be estimated
    New <- 1                                  # subsample (1 = 2008-12 ; 0=2000-04)
 
    vars <- c("pct","POOR_4")
 
    R <- 5000
    matr1 <- data.frame(matrix(0,9,2))

    
    for (i in 1:2){
      for(j in 1:3){
        for(jj in 1:2){
          set.seed(55)
          
          # Set column specifications for Table
          bd <- subset(band.matrix, Var == vars[i] & Kernel == 1)[,1:2] 
          ctr <- 1
          data <- .Clean.Data(M.data,Variable=vars[i],Pov=25,Effect1="State",New=1) 
          formula <- .Formulas(IV=0,controls=1,Variable=vars[i],New=1) 
  
          if(j==1 & jj ==1)  {data <- subset(data, ED.I > 7)                     }
          if(j==1 & jj ==2)  {data <- subset(data, ED.I <= 7)                   }
          if(j==2 & jj ==1)  {data <- subset(data, PA.I %in% parties_nonclient) }
          if(j==2 & jj ==2)  {data <- subset(data, PA.I %in% parties_client)}
          if(j==3 & jj ==1)  {data <- subset(data, PA.I %in% parties_left)}
          if(j==3 & jj ==2)  {data <- subset(data, !(PA.I %in% parties_left))}
  
          ### Cleans Data and Run regressions   
          vec <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                            kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                            reg.end=reg.end, pretre=0,IV=0) 
          output <- .Run.coeff2(vec,reg.beg,reg.end) 
          cit2 <- data.frame(unique(data$City),1) 
          names(cit2) <- c("City","One") 
          
          ### Bootsrap SDs
          results <- boot(data=cit2, statistic=.Boot.coeff, R=R, 
                          base = data, std1 = sdPOP , 
                          std2= sdIDH , kernel = 1, formula = formula,
                          bd = bd, IV = 0, parallel = "snow", 
                          ncpus = 8,reg.end = reg.end, reg.beg = reg.beg, 
                          pop.cut = pop.cut, pretre=0)
          
          output <- .Coeff.intervals(output,results)  
          
          ### prepare output
          OUT <- data.frame(output) 
          
          names(OUT) <- c("COEFF","BCA99L","BCA99H","BCA95L","BCA95H",
                          "BCA90L","BCA90H","Obs","Bpop","Bhdi") 
          OUT$Stars <- ""
          OUT$Stars[(OUT$BCA90L*OUT$BCA90H) > 0] <- "*" ; 
          OUT$Stars[(OUT$BCA95L*OUT$BCA95H) > 0] <- "**" 
          OUT$Stars[(OUT$BCA99L*OUT$BCA99H) > 0] <- "***"
          
         
          matr1[(j-1)*3+1,jj] <- paste(format(round(OUT$COEFF[1], 3),nsmall=3),OUT$Stars[1],sep="")
          matr1[(j-1)*3+2,jj] <- paste("[",format(round(OUT$BCA90L[1], 2),nsmall=2),",",
                                 format(round(OUT$BCA90H[1], 2),nsmall=2),"]",sep="")
          matr1[(j-1)*3+3,jj] <-  paste("{",format(round(OUT$Obs[1], 0),nsmall=0),"}",sep="")
         
        
        }
        
        print(j)
        
      }
      
      if(i==1) {matrx1 <- matr1}else{matrx1 <- cbind(matrx1,matr1)}
   
    }
    
    ccc <-  c("Coefficient","CI","Obs. per bin",
              "Coefficient ","CI ","Obs. per bin ",
              "Coefficient  ","CI  ","Obs. per bin  ")
    table <- .LatexTable(matrx1,ccc,c(3,6,9))
    
    print(table, row.names=FALSE)
    
  #######################################################################################################   
        
  # Figure II ###########################################################################################    
    
 
    docs1 <- dex %>% group_by(Type, name ) %>% summarise_all(c("sum"))
    docs1a <- filter(docs1, Type == "CCT")
    docs1b <- filter(docs1, Type == "Target")
    docs1 <- merge(docs1a,docs1b,by="name")
    docs1 <- docs1 %>% select(name, val.x, val.y)
  
    names(docs1) <- c("name","CCT","Target")

    docs1$Target[nrow(docs1)-2]  <- docs1$Target[nrow(docs1)] 
 
    docs2 <- dex %>% group_by(Type, name ) %>% summarise_all(c("sd"))
    docs2 <- filter(docs2, Type == "Dev")
    docs2$City <- docs2$Type <- NULL
 
    doc <- merge(docs1,docs2,by="name")
    no2 <- doc$CCT[1]/no
    doc$CCT[1:2] <- doc$CCT[1:2] - 974*no2
  
    doc$DEV <- (doc$CCT*1000/doc$Target-1)*100
    doc$name[2] <- "2004H2"
    doc$CCT <- doc$CCT/1000
     
    docx <- filter(doc, name %in% c("Y2004","2004H2","Y2006","Y2008","Y2010","Y2012"))
    docx$Year <- c("2004H1","2004H2","2006","2008","2010","2012")
    
    ggplot(data=docx, aes(x=Year, y=DEV))   + theme_bw()  +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
      geom_bar( stat="identity", position = "dodge", fill="#9ecae1") +
      xlab("Year") + ylab("CCT Coverage, % below target") +
      geom_line(data=docx, aes(x=Year, y=CCT*30-90, group = 1), size=5, color="gray35", alpha=0.3) +
      theme(legend.position="right") +
      theme(axis.text = element_text(colour="black", size=sz)) +
      theme(axis.title = element_text(colour="black", size=sz )) +
      geom_errorbar(aes(ymax = DEV + val, ymin=DEV - val),stat="identity" , width=0.25) +
      geom_hline(yintercept=0)+
      scale_y_continuous("CCT Coverage of poor families, pp",
                         sec.axis = sec_axis(~(.+90)/30, "CCT benefits, mn (gray line)")) 
    
     rm(doc,docs1,docs2,docs1a,docs1b,no,no2)
    
  #######################################################################################################   
  
  # Figure III ##########################################################################################    
 
    bd <- subset(band.matrix, Var == "BF2" & Kernel == 1)[,1:2] 
    data <- .Clean.Data(M.data,Variable="BF2", Pov=25,Effect1="State",New=1) 
    px <- 30 ; py <- 0.7 
    xmin <- 30-1*sdPOP
    xmax <- 30+1*sdPOP
    ymin <- 0.7-1*sdIDH
    ymax <- 0.7+1*sdIDH

    x1<-seq(-1,1,length.out=100) ; y1<- ((1-x1**2)*0.49)**0.5
    x2<-seq(-1,1,length.out=100) ; y2<- -((1-x1**2)*0.49)**0.5
    x <- c(x1,x2) ; y <- c(y1,y2)
    
    sz <- 12 
    test <- data.frame(x=c(seq(from=10,to=30,length=100),seq(from=10,to=30,length=100)),
                       y=c(seq(from=0.65,to=0.75,length=100),seq(from=0.75,to=0.65,length=100)))

    data$D <- as.factor(data$D)
    ggplot(data, aes(x=Pop, y=IDH, color = D)) + theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
      geom_point(shape=1) +
      scale_x_continuous(breaks = round(seq(0, 90, by = 15),1)) +
      theme(axis.text = element_text(colour="black", size=sz) )+
      theme(axis.title = element_text(colour="black", size=sz) )+
      xlab("Population (thousands)") + ylab("HDI") +
      theme(legend.position="none") +
      scale_colour_manual(values=c( "skyblue4","lightskyblue2")) +
      geom_segment(x=2.2,xend=30,y=0.7,yend=0.7,alpha=0.3,size=1,color="red1")+
      geom_segment(x=30,xend=30,y=0.565,yend=0.7,alpha=0.3,size=1,color="red1")
    
  #######################################################################################################   
    
  # Figure IV ###########################################################################################    
 
    reg.beg <- 1
    reg.end <- 19
    bd <- matrix(1,19,2)
    
    for (k in 1:3){
      
      if (k == 1) {Variable <- "ESF"; nom <- "FHP Cov."}
      if (k == 2) {Variable <- "Total_4" ; nom <- "Budget"}
      if (k == 3) {Variable <- "OldCCT"; nom <- "Old CCT Cov."}

      data <- .Clean.Data(M.data,Variable=Variable,Pov=25,Effect1="State",New=1)
      formula <- .Formulas(IV=0,controls=0,Variable=Variable,New=1) 
      if (k == 2) data$Var <- data$Total_4/data$Pop_cur
      vec1 <- .Run.coeff(full.data=data,formula=formula,
                        sdPOP=sdPOP,sdIDH=sdIDH,
                        kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                        reg.end=reg.end, pretre=1, IV=0) 
      

      vec1$avg <- (vec1$X1/sd(data$Var)) 
      vec1$avg <- vec1$avg/vec1$avg[10]
      vec1$rank <- as.factor(c(1:19))
      vec1$name <- nom
  
      
      if(k == 1){jazz <- vec1}else{jazz <- rbind(jazz,vec1)} 
    }

    jazz3 <- jazz[c(1:19,39:57),]
    jazz3$avg <- jazz3$avg+0.6
    
    
    ggplot(data=jazz3 , aes(x=rank, y=avg,  group=name) ) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
      geom_line(data=jazz3 ,aes(x=rank, y=avg, group=name, colour=name), size=2) +
      geom_vline(xintercept = 10, alpha=0.5,size=2,colour="gray75") +
      scale_colour_manual(name="Pre-treatment Averages:",values=c("gray25","gray75")) +
      scale_x_discrete( breaks=c(2,6,10,14,18), 
                        labels=c("Pop=10\n HDI=0.70","Pop=20\n HDI=0.70","Pop=30\n HDI=0.70",
                                 "Pop=30\n HDI=0.65","Pop=30\n HDI=0.60")  )+
      scale_y_continuous(breaks = c(0,0.33,0.67,1)) +
      xlab("") + ylab("") +
      theme(plot.title = element_text(colour="black", size=12)) +
      theme(axis.text = element_text(colour="black", size=10)) +
      theme(legend.text = element_text(colour="black", size = 10)) +
      theme(legend.title = element_text(colour="black", size = 10)) +
      theme(axis.title=element_blank()) +
      theme(legend.position=c(0.25,0.85),legend.direction="vertical") +
      annotate("text",x=3,y=1.7,label="Pre-existing FHP \n coverage is high",size=3.5) +
      annotate("text",x=17,y=1.7,label="Old CCT \n coverage is high",size=3.5) 
    
  #######################################################################################################   
    
  # Figure V ############################################################################################    
    
    reg.beg <- 1
    reg.end <- 19
    
  
    bd <- matrix(1,19,2)
    data <- .Clean.Data(M.data,Variable="BF2", Pov=25,Effect1="State",New=1) 
    formula <- .Formulas(IV=0,controls=1,Variable="BF2",New=1) 
    vec1 <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                       kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                       reg.end=reg.end, pretre=0,IV=0) 
    
    vec1$X10 <- vec1$X10*1000
    jazz <- vec1 ; jazz$rank <- as.factor(c(1:19)) 
    jazz$g <- "a" ; jazz$g[jazz$rank %in% c(9,10,11,12,13,14)] <- "b"
    
    ggplot(data=jazz, aes(x=rank,group=g, y=X1,ymax=max(high),ymin= min(low)) ) + 
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
      geom_vline(xintercept = 10, alpha=0.5,size=2,colour="gray75") +
      geom_point( aes(size = X4, group=g, shape=g),show.legend=FALSE) +
      scale_shape_manual(values=c(1,16)) +
      geom_errorbar(aes(ymin=low, ymax=high), colour="grey15", width=0)+
      scale_x_discrete( breaks=c(2,6,10,14,18), 
                        labels=c("Pop=10\n HDI=0.70","Pop=20\n HDI=0.70","Pop=30\n HDI=0.70",
                                 "Pop=30\n HDI=0.65","Pop=30\n HDI=0.60")  )+
      xlab("") + ylab("") +
      theme(plot.title = element_text(colour="black", size=12)) +
      theme(axis.text = element_text(colour="black", size=10)) +
      geom_hline(yintercept = 0, colour="gray25") +
      theme(axis.title=element_blank()) +
      theme(legend.position="none") 
    
    
  #######################################################################################################   
    
  # Figure VI ###########################################################################################    
    
    reg.beg <- 1
    reg.end <- 19          
    
    nom <- c("lPABV","BF2","pct","Ncand","MV","POOR_4")
    names <- c("Health Funds","CCT Coverage","Incumbent's Share",
               "Number of Candidates","Margin of Victory","Pro-Poor Spending")
    
    
    for (i in 1:6){
      
      bd <- subset(band.matrix, Var == nom[i] & Kernel == 1)[,1:2]
      data <- .Clean.Data(M.data,Variable=nom[i], Pov=25,Effect1="State",New=1) 
      formula <- .Formulas(IV=0,controls=1,Variable=nom[i],New=1) 
      
      vec1 <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                         kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                         reg.end=reg.end, pretre=0,IV=0) 
      vec1$X10 <- vec1$X10*1000
      jazz <- vec1 
      jazz$rank <- as.factor(c(1:19)) 
      
      titl <- names[i]
      
      p <-  ggplot(data=jazz, aes(x=rank, y=X1,ymax=max(high),ymin= min(low)) ) + 
        theme_bw() +
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
        geom_vline(xintercept = 10, alpha=0.5,size=2,colour="gray75") +
        geom_point(shape=16, aes(size = X4), color = "lightskyblue2") +
        geom_errorbar(aes(ymin=low, ymax=high), colour="grey15", width=0)+
        scale_x_discrete( breaks=c(2,6,10,14,18), 
                          labels=c("Pop=10\n HDI=0.70","Pop=20\n HDI=0.70",
                                   "Pop=30\n HDI=0.70",
                                   "Pop=30\n HDI=0.65","Pop=30\n HDI=0.60")  )+
        xlab("") + ylab("") + ggtitle(titl) +
        theme(plot.title = element_text(colour="black", size=12)) +
        theme(axis.text = element_text(colour="black", size=10)) +
        #theme(axis.title.y = element_text(colour="black", size=12, family = "A")) +
        geom_hline(yintercept = 0, colour="gray25") +
        theme(axis.title=element_blank()) +theme(legend.position="none")
      p
      
      if(i == 1) p1 <- p  
      if(i == 2) p2 <- p
      if(i == 3) p3 <- p
      if(i == 4) p4 <- p
      if(i == 5) p5 <- p
      if(i == 6) p6 <- p
      
    }
    
    grid.arrange(p1, p2,p3,p4,p5,p6 ,ncol=2, nrow=3)
    

  #######################################################################################################   
  
  # Figure VII ##########################################################################################    

    reg.beg <- 1
    reg.end <- 19          
    
    names <- c("Challenger's Entry (Education)","Challenger's Success (Education)",
               "Challenger's Entry (Clientelistic)","Challenger's Success (Clientelistic)",
               "Challenger's Entry (Left)","Challenger's Success (Left)")
    nom <- c("chl.entr.educ","chl.succ.educ",
             "chl.entr.client","chl.succ.client",
             "chl.entr.left","chl.succ.left")
    
    for (i in 1:6){
      
      bd <- subset(band.matrix, Var == nom[i] & Kernel == 1)[,1:2]
      data <- .Clean.Data(M.data,Variable=nom[i], Pov=25,Effect1="State",New=1) 
      formula <- .Formulas(IV=0,controls=1,Variable=nom[i],New=1) 
      
      vec1 <- .Run.coeff(full.data=data,formula=formula,sdPOP=sdPOP,sdIDH=sdIDH,
                         kernel=1,bd=bd,pop.cut=pop.cut,reg.beg=reg.beg,
                         reg.end=reg.end, pretre=0,IV=0) 
      vec1$X10 <- vec1$X10*1000
      jazz <- vec1 
      jazz$rank <- as.factor(c(1:19)) 
      
      titl <- names[i]
      
      p <-  ggplot(data=jazz, aes(x=rank, y=X1,ymax=max(high),ymin= min(low)) ) + 
        theme_bw() +
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
        geom_vline(xintercept = 10, alpha=0.5,size=2,colour="gray75") +
        geom_point(shape=16, aes(size = X4), color = "lightskyblue2") +
        geom_errorbar(aes(ymin=low, ymax=high), colour="grey15", width=0)+
        scale_x_discrete( breaks=c(2,6,10,14,18), 
                          labels=c("Pop=10\n HDI=0.70","Pop=20\n HDI=0.70",
                                   "Pop=30\n HDI=0.70",
                                   "Pop=30\n HDI=0.65","Pop=30\n HDI=0.60")  )+
        xlab("") + ylab("") + ggtitle(titl) +
        theme(plot.title = element_text(colour="black", size=12)) +
        theme(axis.text = element_text(colour="black", size=10)) +
        #theme(axis.title.y = element_text(colour="black", size=12, family = "A")) +
        geom_hline(yintercept = 0, colour="gray25") +
        theme(axis.title=element_blank()) +theme(legend.position="none")
      p
      
      if(i == 1) p1 <- p  
      if(i == 2) p2 <- p
      if(i == 3) p3 <- p
      if(i == 4) p4 <- p
      if(i == 5) p5 <- p
      if(i == 6) p6 <- p
      
    }
    
    grid.arrange(p1, p2,p3,p4,p5,p6 ,ncol=2, nrow=3)
    

  #######################################################################################################   
    
    
